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Abstract 

We investigate the controversial issue of the existence of universality classes 
describing critical phenomena in three-dimensional systems characterized by a 
matrix order parameter with symmetry 0(2)<8>0(iV) and symmetry-breaking 
pattern 0(2)(g>0(iV)^0(2)(5<)0(iV — 2). Physical realizations of these systems 
are, for example, frustrated spin models with noncollinear order. 

Starting from the field-theoretical Landau-Ginzburg- Wilson Hamiltonian, 
we consider the massless critical theory and the minimal-subtraction scheme 
without e expansion. The three-dimensional analysis of the corresponding 
five-loop series shows the existence of a stable fixed point for N = 2 and N = 
3, confirming recent field-theoretical results based on a six-loop expansion in 
the alternative zero-momentum renormalization scheme defined in the massive 
disordered phase. 

In addition, we report numerical Monte Carlo simulations of a class of 
three-dimensional 0(2)(g>0(2)-symmetric lattice models. The results provide 
further support to the existence of the 0(2)<g)0(2) universality class predicted 
by the field-theoretical analyses. 

PACS Numbers: 05.70.Jk, 75.10.-b, 05.10.Cc, 64.60.Fr 
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I. INTRODUCTION 



Several interesting critical transitions are effectively described by a matrix order parame- 
ter with symmetry 0(2)®0(iV) and symmetry-breaking pattern 0(2)®0(iV)-^0(2)(g>0(iV— 
2). This is the case, for N = 2 or N = 3, of multicomponent frustrated magnetic systems 
with noncollinear order, in which frustration may arise either because of the special geome- 
try of the lattice or from the competition of different kinds of interactions. Typical examples 
of systems of the first type are stacked triangular antiferromagnets (STA's), in which the 
magnetic ions are located at the sites of a stacked triangular lattice. Frustration due to the 
competition of interactions is realized in helimagnets, in which a magnetic spiral is formed 
along a certain direction of the lattice. The nature of the magnetic transition in these ma- 
terials has been the object of several studies, see, e.g., Refs. [1-5] for reviews. In particular, 
the order of the transition is still controversial, with several contradictory results both on 
theoretical and experimental sides. 

The Landau-Ginzburg- Wilson (LGW) theory with 0(2)®0(iV) symmetry that is ex- 
pected to describe these systems is given by [6] 

^LGW = J d d x{^J2[( d M 2 + r<t>^ +^o(^^) 2 + ^oJ][(0a-0 6 ) 2 -0^]}, 

"a a a,b 

(1.1) 

where <p a {a = 1, 2) are iV-component vectors. The symmetry-breaking pattern [7] 
0(2)<g>0(iV)— >0(2)(g)0(./V — 2) is obtained by requiring vq > 0. Negative values of vq lead 
to a different symmetry-breaking pattern: the ground state configurations have a ferromag- 
netic or antiferromagnetic order and correspond to 0(2)®0(A r )^ Z 2 ®0(iV — 1). For v < 
the model is also of interest; it describes magnets with sinusoidal spin structures [8,6] and, 
for N = 3, the superfluid transition of 3 He [9-11]; see, e.g., Refs. [12,3,4,13] for other ap- 
plications. Here, we will only focus on the case v o > and thus whenever we speak of an 
0(2)(g)0(iV) universality class we refer to the case in which the symmetry-breaking pattern 
is 0(2)<g>0(iV)->0(2)<8>0(iV - 2). 

The 0(2)£g>0(iV) theory (1.1) has been much studied using field-theoretical (FT) meth- 
ods. Different perturbative schemes have been exploited, such as the e = 4 — d expansion [14] 
and the three-dimensional (3-d) massive zero-momentum (MZM) renormalization scheme 
[15]. A detailed discussion of the scenario emerging from the e expansion is presented in 
Ref. [3]. Near four dimensions a stable 0(2)®0(iV) fixed point (FP) with v > is found 
[6,16,17] only for large values of N, N > N c = 21.80-23.43e + 7.09e 2 + O(e 3 ). Resummations 
of the e expansion of N c , known to 0(e 4 ) [17], suggest [17,18] N c m 6 in three dimensions. 
Therefore, a smooth extrapolation of the scenario around d = 4 to d = 3 would indicate 
that a new 0(2)®0(iV) universality class does not exist for the physically interesting cases 
N = 2 and 3. On the other hand, six-loop calculations in the framework of the 3-d MZM 
scheme provide a rather robust evidence for the existence of a new stable FP for iV = 2 
and N = 3 with attraction domain in the region v o > [19,20]. This FP was found only 
in the analysis of high-order series, starting at four loops, while earlier lower-order calcu- 
lations up to three loops [12] did not find it. According to renormalization-group (RG) 
theory, the stable FP of the 0(2)®0(iV) theory should describe the critical behavior of 3-d 
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systems undergoing continuous transitions characterized by the symmetry-breaking pattern 
0(2)®0(iV)->0(2)<8>0(iV - 2). The main problem of the calculations within the MZM 
scheme is the fact that, for N = 2 and 3, the 0(2)(g>0(iV) FP is found in a region of quartic 
couplings in which the perturbative expansions are not Borel summable. Therefore, a Borel 
transformation only provides an asymptotic expansion and convergence is not guaranteed, at 
variance with the case of O(N) theories in which the Borel summability of the corresponding 
MZM expansion provides a solid theoretical basis for the resummation methods. In the case 
of the 0(2)<g)0(iV) theory the reliability of the results concerning the new stable FP is es- 
sentially verified a posteriori from their stability with respect to the perturbative order. The 
MZM expansions have also been analyzed by using the pseudo-e expansion method [17,21]. 
No stable FP is found for N = 2 and 3, but this is not unexpected since this resummation 
method can only find FPs that are already present at one loop, similarly to the e expansion. 
We finally mention that perturbative studies of the corresponding nonlinear a models near 
two dimensions have been reported in Refs. [22,23]. 

We mention that there are other physically interesting cases in which low-order e- 
expansion calculations fail to provide the correct physical picture: for example, the Ginzburg- 
Landau model of superconductors, in which a complex scalar field couples to a gauge field. 
Although e-expansion calculations do not find a stable FP [24], thus predicting first-order 
transitions, it is now well established (see, e.g., Refs. [25,26]) that 3-d systems described 
by the Ginzburg-Landau model can also undergo a continuous transition — this implies the 
presence of a stable FP in the 3-d Ginzburg-Landau theory — in agreement with experiments 
[27]. 

The 0(2)(g)0(iV) theory has also been studied by exploiting an alternative FT method 
based on the analysis of the RG flow of the so-called effective average action [28-31,5]. 
This approach does not rely on a perturbative expansion around the Gaussian FP and it is 
therefore intrinsically nonperturbative. However, the practical implementation requires ap- 
proximations and truncations of the average effective action. For this purpose, a derivative 
expansion of the effective average action is usually performed. The studies of the 0(2)®0(iV) 
theory reported in the literature [29-31,5], based on the zeroth- and first-order approxima- 
tions, do not find evidence of stable 0(2)®0(iV) FPs for iV = 2 and 3, in contradiction with 
the perturbative MZM results. This would imply that phase transitions characterized by 
the symmetry-breaking pattern 0(2)®0(iV)— >0(2)®0(iV — 2) with N = 2 or 3 are always 
of first order. 

The issue concerning the existence of 0(2)®0(iV) universality classes is most important 
to understand the physics of STA's and of magnets with helical order, because the absence of 
stable 0(2)(g)0(iV) FPs implies that none of them can undergo a continuous transition. On 
the experimental side, experiments [1,2] have apparently observed continuous transitions be- 
longing to the 0(2)cg>0(iV) universality class. However, as discussed in Ref. [5], experimental 
results are not consistent — STA's and helimagnets show a critical behavior with apparently 
different exponents — and, in some cases, do not satisfy general exponent inequalities, for 
instance 7 < 2v and 2(3 > v. 

The most recent Monte Carlo (MC) simulations of the antiferromagnetic XY Hamil- 
tonian on a stacked triangular lattice have observed a first-order transition [32-34] with 
very small latent heat. Moreover, first-order transitions have been observed in MC investi- 
gations [35,34] of modified lattice spin systems whose transitions are characterized by the 
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same symmetry-breaking pattern. Therefore, MC simulations of the models considered up 
to now do not support the existence of an 0(2)<g>0(2) universality class. On the other hand, 
MC simulations of Heisenberg STA models, corresponding to N = 3, give results that are 
substantially consistent with a continuous transition, see, e.g., Ref. [36]. 

We would like to stress that the existence of a universality class is not contradicted by the 
observation of first-order transitions in some systems sharing the same symmetry-breaking 
pattern. The universality class determines the critical behavior only if the system undergoes 
a continuous transition. Instead, first-order transitions are expected for systems that are 
outside the attraction domain of the stable FP. This is evident in mean-field calculations and 
also within the FT approach, in which some RG trajectories do not flow towards the stable 
FP but run away to infinity. Therefore, the above-mentioned MC results for the XY STA 
models are still compatible with the hypothesis of the existence of an 0(2) (80 (2) universality 
class; XY STA models may be simply outside the attraction domain of the stable FP. 

In this paper we further investigate the existence of the 0(2)<g>0(iV) universality class 
for XY (N = 2) and Heisenberg (N = 3) systems. First, we consider an alternative 3-d 
perturbative approach, the so-called minimal-subtraction (MS) scheme without e expansion 
[37-39], for which five-loop series have been recently computed in Ref. [17]. This scheme is 
strictly related to the e expansion, but, unlike it, no e expansion is performed and e is set 
to the physical value e = 1, providing a 3-d scheme. It works within the massless critical 
theory, thus providing a nontrivial check of the results obtained within the MZM scheme, 
which is defined in the massive disordered phase. The analysis of the corresponding five- 
loop expansions shows the existence of an 0(2)(g)0(iV) FP for N = 2 and 3, confirming the 
conclusions of the analysis of the six-loop expansions within the MZM scheme. Concerning 
the critical exponents, the analysis of the five-loop MS series gives v = 0.65(6) and r) = 
0.09(4) for N = 2, and v = 0.63(5) and rj = 0.08(3) for N = 3. These results should be 
compared with the estimates obtained from the six-loop MZM series, which are v = 0.57(3) 
and 7] = 0.09(1) for N = 2, and v = 0.55(3) and rj = 0.10(1) for N = 3. It is important to 
note that, although the available MS series have one order less, the corresponding results are 
expected to be more reliable than the MZM ones, because the MS FPs are at the boundary 
of the region in which the expansions are Borel summable, and not outside it as in the 
MZM case. We finally mention that the MS scheme without e expansion allows us to obtain 
fixed-dimension results at any d. Thus, we are able to recover the results of the e expansion 
sufficiently close to four dimensions and to obtain a full picture of the fate of the different 
FPs as d varies from four to three dimensions. 

We also address numerically the question of identifying a 3-d lattice model with symmetry 
0(2) <8>0 (2) and with the expected symmetry-breaking pattern that shows a continuous 
transition. This would conclusively show that the 0(2) (80 (2) universality class really exists. 
For this purpose we consider the following lattice model 

w = -/3 Y fax ■ ¥x+n + ipx ■ ifrx+n) + Y Wx + + 

X,fl X 

+A, £ - If + (^ - I) 2 ] + 2A 22 &I>1 (1-2) 

X X 

where tp and ip are two-component real variables. The Hamiltonian 7i describes two identical 
two-component 0(2)-symmetric lattice 4 models coupled by an energy-energy term. By an 
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appropriate change of variables, see Sec. IV, one can show that model (1.2) corresponds to a 
lattice discretization of the Hamiltonian (1.1) for N = 2 with uq ~ (^4+^22) and Vq ~ (A 2 2 — 
A 4 ). When A 22 > A 4 > the critical behavior at the phase transition should be described 
by the N = 2 Hamiltonian (1.1) with v > 0. Therefore, a region of continuous transitions 
in the quartic parameter space with A22 > A4 would imply the existence of the 0(2)<g>0(2) 
universality class. In order to investigate this point, we present MC simulations for A4 — I 
and several values of ^22- The phase diagram emerging from the simulations is characterized 
by a line of first-order transitions extending from large values of A22 down to a tricritical 
point at A 2 2 = A* 22 > A 4 — 1, where the latent heat vanishes, and, for 1 = A 4 < A 2 2 < A 22 , 
by a line of continuous transitions that should belong to the 0(2)®0(2) universality class 
identified by the perturbative FT approaches. The possible extension of the first-order 
transition line up to A22 = A4, i.e. up to the 4-vector theory, is apparently incompatible 
with the theoretically predicted behavior of the latent heat near an 0(4) tricritical point. 

The paper is organized as follows. In Sec. II we present the analysis of the five-loop MS 
expansions, providing evidence for the existence of a stable FP with attraction domain in the 
region vq > 0, in the two- and three-component cases. There, we also show that for d — > 4 
the results of the e expansion are recovered. In Sec. Ill we discuss the crossover behaviors 
predicted by the FT approach and their relation with those that may be observed in realistic 
models. In Sec. IV we report the results of the MC simulations for the model defined 
by Hamiltonian (1.2), determining its phase diagram in the region of quartic parameters 
A22 > At = 1- We investigate the finite-size scaling (FSS) behavior using cubic lattices 
of size L < 120. In Sec. V we report some conclusive remarks. In App. A we provide 
some details on the perturbative expansions in the MS scheme. In App. B we discuss some 
properties of 0(M)®0(V)-symmetric medium-range models. 

II. ANALYSIS OF THE FIVE-LOOP MINIMAL-SUBTRACTION SERIES 

A. The minimal-subtraction scheme without e expansion 

The FT approach is based on the Hamiltonian (1.1). In the MS scheme one considers the 
massless critical theory in dimensional regularization and determines the RG functions from 
the divergences appearing in the perturbative expansion of the correlation functions [37]. In 
the standard e-expansion scheme [14] the FPs, i.e. the common zeroes of the /3-functions, 
are determined perturbatively as expansions in powers of e, while exponents are obtained by 
expanding the corresponding RG functions, i.e. 77^ (see App. A), computed at the FP in 
powers of e. The MS scheme without e expansion [39] is strictly related. The RG functions 
P u>v and rj^t are the MS functions. However, e is no longer considered as a small quantity 
but it is set to its physical value, i.e. in three dimensions one simply sets e = 1. Then, 
one determines the FP values u*, v* from the common zeroes of the resummed f3 functions. 
Finally, critical exponents are determined by evaluating the resummed RG functions r]^ and 
f) t at u* and v*. Notice that the FP values u* and v* are different from the FP values of the 
renormalized quartic couplings of the MZM renormalization scheme, since u and v indicate 
different quantities in the two schemes. 
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The MS RG functions have been computed to five loops in Ref. [17]. In App. A we report 
the series for N = 2 and N = 3. We also consider the critical exponents associated with 
the chiral degrees of freedom. They can be determined from the RG dimension of the chiral 
operator [6] 

C ck di(x) = <frck(x)(j)di(x) - 4>ci{x)(t>dk{x). (2.1) 

We computed the MS RG function rj c (u : v) associated with the chiral operator C c kdi to five 
loops. The series are reported in App. A. 



B. The resummation 

Since perturbative expansions are divergent, resummation methods must be used to 
obtain meaningful results. Given a generic quantity S(u, v) with perturbative expansion 
S(u,v) = CijU % v^ , we consider 

S(xu, xv) = ^2 s k(u, v)x k , (2.2) 
k 

which must be evaluated at x — 1. The expansion (2.2) in powers of x is resummed by 
using the conformal-mapping method [40,41] that exploits the knowledge of the large-order 
behavior of the coefficients, generically given by 

s k (u, v) ~ k\ [-A(u, v)] k k b [1 + Oik' 1 )] . (2.3) 

The quantity A(u,v) is related to the singularity t s of the Borel transform B(t) that is 
nearest to the origin: t s = —\/A{u,v). The series is Borel summable for x > if B(t) does 
not have singularities on the positive real axis, and, in particular, if A(u,v) > 0. Using 
semiclassical arguments, one can argue that [19] the expansion is Borel summable when (see 
App. A for the precise definition of u and v) 

m>0, u-^v>0. (2.4) 

In this region we have 

A(u,v) = ^Max(u,u-v/2). (2.5) 

Under the additional assumption that the Borel-transform singularities lie only in the nega- 
tive axis, the conformal-mapping method [40] turns the original expansion into a convergent 
one in the region (2.4). Outside, the expansion is not Borel summable. However, if the 
condition 

u-^v>0 (2.6) 

holds, then the Borel-transform singularity closest to the origin is still in the negative axis, 
and therefore the large-order behavior is still given by Eq. (2.3) with A(u,v) given by 
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Eq. (2.5). Thus, by using this value of A(u,v) and the conformal-mapping method one 
still takes into account the leading large-order behavior. One may therefore hope to get an 
asymptotic expansion with a milder behavior, which may still provide reliable results. 

We should mention that the MS series are essentially four-dimensional, so that they may 
be affected by renormalons that make the expansion non-Borel summable for any u and 
v, and are not detected by a semiclassical analysis; see, e.g., Ref. [42]. The same problem 
should also affect the MS series of O(N) models. However, the good agreement between the 
results obtained from the FT analyses [39] and those obtained by other methods indicates 
that renormalon effects are either very small or absent (note that, as shown in Ref. [43], 
this may occur in some renormalization schemes). For example, the analysis of the five-loop 
perturbative series [39] gives v = 0.629(5) for the Ising model and v = 0.667(5) for the 
XY model, that are in good agreement with the most precise estimates obtained by lattice 
techniques, such as v = 0.63012(16) [44] and v = 0.63020(12) [45] for the Ising model, and 
v = 0.67155(27) [46] for the XY universality class. On the basis of these results, we will 
assume renormalon effects to be negligible in the analysis of the two-variable series of the 
0(2) (8)0 (iV) theory. 

C. Three-dimensional analysis of the five-loop series for N = 2,3 

The RG flow of the theory is determined by the FPs. Two FPs are easily identified: the 
Gaussian FP, which is always unstable, and the 0(2iV) FP located along the u axis. The 
results of Ref. [48,49] on the stability of the three-dimensional 0(M)-symmetric FP under 
generic perturbations can be used to prove that also the 0(2N) FP is unstable for any N > 2 
[47]. Indeed, the Hamiltonian term (0 a • <f)b) 2 — which acts as a perturbation at the 

0(2N) FP, is a particular combination of quartic operators transforming as the spin-0 and 
spin-4 representations of the 0(2iV) group, and any spin-4 quartic perturbation is relevant 
[48] at the O(M) FP for M > 3, since its RG dimension y AA is positive for M > 3. In 
particular, y 4A « 0.11 at the 0(4) FP and y AA « 0.27 at the 0(6) FP [48]. Note that these 
values are rather small, especially in the 0(4) case. The w-axis plays the role of a separatrix 
and thus the RG flow corresponding to v > cannot cross the u-axis. Therefore, since 
models with the symmetry-breaking pattern 0(2)®0(iV)-^0(2)£g>0(iV — 2) have t>o > 0, the 
relevant FPs lie in the region v > 0. 

The analyses of the six-loop series in the MZM scheme reported in Refs. [19,20] provided 
a rather robust evidence for the presence of a stable FP with attraction domain in the region 
v o > for N = 2 and N = 3. In the following this result will be confirmed by the analysis 
of the MS five-loop series. In order to investigate the RG flow in the region v > 0, we 
apply essentially the same analysis method of Refs. [19,50] (we refer to these references for 
details). We resum the perturbative series by means of the conformal-mapping method and, 
in order to understand the systematic errors, we vary two different parameters b and a (see 
Ref. [50] for definitions). We also apply this method for those values of u and v for which 
the series are not Borel summable but still satisfy u — \v > 0. As already discussed, the 
conformal-mapping method should still provide reasonable estimates since we are taking 
into account the leading large-order behavior. 

In order to find the zeroes of the /^-functions we first resummed the expansions of B u (u, v) 
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FIG. 1. Zeroes of the (5- functions in the region v > for N = 2 and 3 for a particular 
approximant, see text. The three colored regions correspond (from below) to: (1) u — \v > 
(domain in which the perturbative series are Borel summable), (2) u — \v > and u — \v < 
(domain in which the perturbative series are not Borel summable but one can take into account 
the leading large-order behavior), (3) u — \v < 0. The full lines at the 0(2)(g>0(iV) FP show the 
final estimate (2.7) with its uncertainty. 



and B v (u, v) denned in Eq. (A4). More precisely, we considered the functions R u ^ v (u, v, x) = 
B U}V (ux, vx)/x 2 . For each function R UyV we considered several approximants corresponding 
to different values of the resummation parameters a and b, see, e.g., Refs. [50,19] for details. 
In Fig. 1 we show the common zeroes of the /3-functions in the region v > 0. The figure is 
obtained by using a single approximant, the one with a u = a v = 1, b u = b v = 10, but others 
give qualitatively similar results. A common zero (u*,v*) with v * > is clearly observed at 
u* w 1.1, v* w 2.3 for N = 2, and at u* ps 0.9, v* w 1.7 for = 3. In order to give an 
estimate of the FP, we considered resummations of B u (u,v) and B v (u,v) with parameters 
a u , b u , a v , and b v , assuming integer values in the range — 1 < a UyV < 3 and 4 < b UyV < 16. 
Most combinations, approximately 90% for N = 2 and 97% for N = 3, have a common 
zero in the region v > (these percentages increase if we only consider approximants with 
ot u = ot v and b u = b v , becoming approximately 94% for N = 2 and 99% for N = 3). We 
take the average of the results as final estimate, obtaining 

u* = 1.10(13), v* = 2.30(21) for N = 2, 

u* = 0.90(8), v* = 1.72(15) for N = 3. (2.7) 

The errors are related to the variation of the results with respect to changes of the resumma- 
tion parameters a u , b u , a v , and b v in the considered range of values, and correspond to one 
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standard deviation. As a check, we also tried a different method. We determined optimal 
values of a and b by minimizing the difference between the results of the four- and five-loop 
resummations of the functions B u v (independently) close to the 0(2)(g)0(iV) FP. The re- 
sults are consistent with those reported in Eq. (2.7). Notice that in the case N = 3, since 
v*/u* ~ 1.9, the FP is substantially within the region in which the perturbative expansions 
should be Borel summable, while for N = 2, since v*/u* ~ 2.1, the FP is slightly above 
its boundary v/u = 2. Therefore, Borel resummations are expected to be effective. In this 
respect the MS scheme seems to behave better than the MZM scheme, in which the FPs are 
in the non-Borel summable region [19], although still in the region in which the conformal- 
mapping resummation method should be able to take into account the leading large-order 
behavior. The analysis of the stability matrix shows that the FP is stable, i.e. its eigen- 
values have positive real part. Most approximants give complex eigenvalues, supporting the 
hypothesis that the FP is a focus, as discussed in Ref. [20]. We obtain u = 1.0(5) ± i0.8(5) 
for N = 2 and uo = 0.9(4) ± i0.7(3) for N = 3, in rough agreement with the MZM scheme 
results [20]. We finally mention that consistent results are obtained by resumming the series 
using the Pade-Borel technique, which does not exploit the knowledge of the large-order 
behavior of the series. 

Critical exponents are obtained by evaluating the RG functions rj^ and r] t or appropriate 
combinations at the FP. We found 

v = 0.65(2 + 4), r] = 0.09(2 + 2), 7 = 1.24(3 + 8), <p c = 1.42(6 + 10) (2.8) 

for N = 2, and 

v = 0.63(1 + 4), 77 = 0.08(2 + 1), 7 = 1.20(1 + 7), <p c = 1.35(2 + 7) (2.9) 

for N = 3. The errors are reported as the sum of two terms, related respectively to the 
dependence on b and a and to the uncertainty of the FP coordinates. For comparison we 
report the corresponding results obtained from the analysis of the six-loop series in the 
MZM scheme [19,51]: v = 0.57(3), rj = 0.09(1) and <p c = 1.43(4) for N = 2, and v = 0.55(3), 
rj = 0.10(1) and C = 1.27(4) for N = 3. We note that the MS estimates of v and 7 
are larger than those obtained in the MZM scheme, but still substantially compatible with 
them taking into account their relatively large errors. We stress again that this comparison 
represents a nontrivial consistency check since the two schemes are quite different: in the 
MZM scheme one works in the massive high-temperature phase, while in the MS scheme 
one considers the massless critical theory. 

Finally, we computed the RG flow, in order to determine the attraction domain of the 
stable FP. We refer the reader to Ref. [52] for the relevant definitions. In Fig. 2 we show 
the RG flow in the quartic-coupling u, v plane corresponding to different values of the ratio 
s = Vq/uq for v > 0. All trajectories corresponding to s < 3/2 belong to the region 
u — \v > 0, in which the resummation should be reliable, and appear to be attracted by the 
stable FP. 

D. Results for generic values of N and dimension d 

The fixed- dimension MS scheme allows us to obtain fixed-dimension results for any di- 
mension d. Since in three dimensions and for = 2, 3 this scheme provides results that 
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u 

FIG. 2. The RG flow in the quartic-coupling plane for N = 2, 3 in the MS scheme for d = 3. 

are substantially different from those of the strictly related e expansion, it is interesting to 
compare the two perturbative methods for generic values of d and N. 

Using the five-loop series of the MS /3-functions, we investigate the presence of FPs in 
the region < v < 4u, where resummations seem to be under control, for generic d and N. 
In any d, or e = 4 — d, and for sufficiently large values of N, we find a stable 0(2)(8>0(iV) 
FP. If we decrease N at fixed e, for e smaller than a critical value e Cjmax , we find a value 
iV c (e) such that, for N = N c (e), the stable FP disappears. In Fig. 3 we plot the results for 
the inverse quantity e c (N) = 4 — d c (N), where d c (N) represents the dimension below which 
one finds a stable FP in the region v > at fixed N. This quantity may be estimated by 
averaging the values of the largest dimension d (smallest value of e) for which each pair of 
approximants of the /3-functions (we use the same set as in the three-dimensional analysis 
reported in Sec. II C) has a stable 0(2)®0(iV) FP. The reported error corresponds to one 
standard deviation. Unfortunately, for 4 < N < 5 the stable FP moves outside the region 
u — \v > 0, in which we are able to resum reliably the perturbative series (if this condition 
is satisfied we can take into account the leading large-order behavior). Therefore, for N < 4 
we are unable to compute e c (N). In this case we can compute a conservative upper bound 
by finding the smallest value of e such that at least 95% of the approximants still present 
a stable FP in the region < v < Au. The bounds corresponding to N = 2, 3, and 4 
are represented by thick segments in Fig. 3. There, we also compare these results with 
the curve obtained by resumming the 0(e 4 ) expansion [17] of N c (e) (we actually report the 
curve obtained by resumming 1/N c (e)). A nice agreement is observed for sufficiently large 
values of N, down to iV w 8. For smaller values of N the fixed-dimension MS results differ 
from the e-expansion curve. In particular, unlike the e expansion, the fixed-dimension MS 
series provides estimate for e c (2) and e c (3) that are definitely smaller than one, leading to the 
bounds e c (2) < 0.7 and e c (3) < 0.8. As Fig. 3 shows, the results for e c (iV) are nonmonotonic, 
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FIG. 3. Results for e c (N) as obtained from the fixed-dimension MS analysis and from the e 
expansion (dashed line). The thick segments at N = 2,3,4 represent conservative upper bounds 
on e c (N). 



with a maximum value e Cjmax ~ 0.8 for N rs 8. Thus, for e < e Cjmax there exists another 
limiting value of N, N c2 (e), such that for N c2 (e) < N < iV c (e) no stable FP exists in the 
region v > 0, while for N < N c2 (e) the 0(2)®0(iV) FP is again present. Note that while for 
iV > N c (e) the stable FP has real stability eigenvalues, for N < N c2 (e) it is a stable focus, 
i.e., the stability eigenvalues are complex with positive real part. 

The MS results are qualitatively consistent with the MZM ones, although the estimate 
of e c ,max, e c ,max ~ 0.8, apparently contradicts the conclusions of Ref. [19]. Indeed, from the 
analysis of the MZM scheme expansions, Ref. [19] did not find a clear evidence of stable 
FPs for 5 < iV < 7, which would imply that e Cjmax > 1. This conclusion was also in contrast 
with the MC results of Ref. [53] that apparently found a continuous transition for N = 6. 
On the basis of the present analysis we are now inclined to believe that this may be only 
a resummation problem, in some sense connected to the fact that in the extended space 
(N, d), N = 6 and d — 3 is close to the line d = d c (N). 

Although in the three-dimensional MS calculation, N = 6 does not represent a special 
point for the existence of the FP, such a value still plays a special role. In d — 3 the FP is 
topologically different for small and large values of N since the stability eigenvalues are real 
for large N, while for N = 2 and 3 they are complex. Therefore, there should be a value 
iVeq that separates the two behaviors. We shall show that N cq w 6. 

For large values of iV the stability eigenvalues are real. As iV decreases the difference 
between u>i and ui 2 decreases and for iV = iV eq we have uj\ = u 2 . Then, for iV < iV eq the 
eigenvalues become complex and the FP is a focus. As it can be seen from the results 
reported in Table I, N eq rs 6. In this case, ~ 50% of the approximants give real estimates 
for ui and u 2 , while ~ 50% give complex estimates with a small imaginary part. In all cases 
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TABLE I. Results obtained from the analysis of the five-loop series in the MS perturbative 
scheme. No values of uj± t 2 are reported for N = 6: in this case it is not clear whether the eigenvalues 
of the stability matrix are complex or real. 



N 


u* 


V* 




V 


7 


V 


2 


1.10(13) 


2.30(21) 


1.0(5) ±i0.8(5) 


0.65(6) 


1.24(11) 


0.09(4) 


3 


0.90(8) 


1.72(15) 


0.9(4) ±i0.7(3) 


0.63(5) 


1.20(8) 


0.08(3) 


4 


0.74(3) 


1.29(8) 


0.7(2) ±i0.4(2) 


0.64(4) 


1.24(6) 


0.073(10) 


5 


0.63(3) 


1.04(7) 


0.7(2) ±i0.3(2) 


0.64(4) 


1.25(6) 


0.061(10) 


6 


0.56(4) 


0.86(7) 




0.66(4) 


1.29(8) 


0.052(14) 


7 


0.51(5) 


0.73(5) 


0.8(2) 0.5(2) 


0.68(4) 


1.34(9) 


0.047(15) 


8 


0.47(4) 


0.64(4) 


0.8(2) 0.5(2) 


0.70(5) 


1.37(10) 


0.042(10) 


10 


0.41(4) 


0.51(1) 


0.9(2) 0.5(2) 


0.74(5) 


1.46(11) 


0.036(8) 


16 


0.300(14) 


0.334(8) 


0.9(2) 0.74(12) 


0.82(4) 


1.62(8) 


0.025(4) 



the real part satisfies 0.3 < StecUj < 0.8. 

Beside the stability eigenvalues in Table I we also report the critical exponents and the FP 
coordinates for several values of N. These results are in agreement with the MZM estimates 
of Ref. [20]. We also note that the MS results for N = 6 are in substantial agreemeent with 
the MC results of Ref. [53], v = 0.700(11) and 7 = 1.383(36), and with the nonperturbative 
RG results of Ref. [5], v « 0.707 and 7 « 1.377. 



III. CROSSOVER BEHAVIOR 
A. Effective exponents 

The perturbative analysis presented in Sec. II as well as the analyses in the MZM scheme 
of Ref. [19] predict the presence of a stable FP for the physically interesting cases N = 2 
and N = 3. However, this FP has a quite unusual feature: the stability eigenvalues are 
apparenly complex with positive real part [19,20]. In this Section we wish to understand the 
consequences on experimental and numerical determinations of the critical exponents and, 
in general, of RG-invariant quantities. 

The presence of complex stability eigenvalues changes the approach to criticality. If O 
is a generic critical quantity we expect close to the critical point 

O « CC [l + a£- w *cos(c;/logf + &)] , (3.1) 

where £ is the correlation length and the stability eigenvalues are written as ujR^siiujj. Scaling 
corrections oscillate and the approach to the asymptotic behavior is nonmonotonic. 

In order to characterize the behavior of critical quantities outside the critical point, it is 
useful to introduce effective exponents. From the susceptibility x an d the correlation length 
£ one can define the effective exponents 

u\ _ 91n £ <91nx f . d\nx f . 

where t = (T — T c )/T c is the reduced temperature. One can easily check that 7] c q = 
2 — 7eff/ i/ eff- The effective exponents are not universal and depend on the specific model. 
Nonetheless, it is usually assumed (but there are notable exceptions; for instance, the 3-d 
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Ising model and the corresponding scalar 4 theory behave differently near the critical point 
[54-56]) that the qualitative features are similar in all models belonging to the same uni- 
versality class. For this reason, in the following we shall compute the effective exponents in 
the FT model. We shall present numerical results for N = 2, in order to be able to compare 
them with the MC results of Sec. IV. For N = 3 effective exponents are qualitatively similar. 

We shall consider the MZM scheme since all necessary formulas have already been pre- 
sented in Ref. [52], although the same analysis could have been done in the MS scheme by 
generalizing to the present case the results of Ref. [57]. If u and v are the zero-momentum 
renormalized couplings normalized so that u ~ Uq/tji and v ~ Vq/tu at tree level [58], RG 
trajectories are determined by solving the differential equations 

nil 
dX 

on 

-A- = /WA), W (A)), (3.3) 

where A G [0, oo), with the initial conditions 

u(0) =v(0) = 0, 
du 



dX 



A=0 



dv 
dX 



(3.4) 



A=0 



where s = Vq/uq parametrizes the different models. The results of Ref. [52] allow us to 
derive general scaling formulas for the rescaled x = X u o an d £ = £mo> where x an d £ are 
respectively the susceptibility and the second-moment correlation length. In particular, if 
r = r — r c is the reduced temperature and r = r/wjj, we have 
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X = F x (t,s), Z = F c (t,s). (3.5) 

The functions F x and can be expressed in terms of RG functions — in the present case they 
are known to six loops — and can be computed rather accurately, as we shall show below. 

In Fig. 4 we show the RG trajectories for several values of s with < s < 1 [59]. For 
larger values of s trajectories run in the region u — \v < 0, where we are not able to resum 
the perturbative series. Correspondingly, in Fig. 5 we report the behavior of the four-point 
couplings u and v as a function of £. The corresponding FP values [19,58] are u* = 23.9(1.3) 
and v* = 68.7(2.5). Considering first v(£), it is interesting to observe that oscillations are 
significant only for s ~ 1. For smaller values of s, v(£) increases essentially monotonically 
with £. More peculiar is the behavior of u(£). Indeed, for all s < 2/3, u flattens first at a 
value around 18 and then suddenly increases towards the asymptotic value. This is due to 
the presence of the unstable 0(4) FP that gives rise to strong crossover effects, even when 
s is as large as 2/3. Indeed, the plateau observed in u corresponds to the FP value of u in 
the 0(4) theory [60], u* Q ^ ~ 17.4. Thus, unless s ~> 1, the flow first feels the presence of 
the 0(4) FP, so that u « «o(4) anc ^ then S oes towards the 0(2)®0(2) FP. In Fig. 6 we also 
show the ratio v/u. Note that for small s such a ratio is very small, while for s > 1/2 the 
behavior is nonmonotonic with a pronounced peak. 

Finally, we determine the effective exponents. In Fig. 7 we report the effective exponents 
u cS and 7/efr as a function of the rescaled reduced temperature f. The exponent z/ eff shows 
quite large oscillations, especially for small s. They are not only due to the complex stability 
eigenvalues but also to crossover effects related to the presence of the O (4) FP. As we already 
remarked above, for small s the trajectories are close to the 0(4) FP and thus u c g is close 
to z/o(4) (the best available estimate is vo{4) = 0.749(2), Ref. [61]). For instance, for s = 1/4 
(resp. s = 1/2) the maximum value of v e s is 0.71 (resp. 0.67). For s close to 1, crossover 
effects are less relevant and z/ e g does not increase much. However, in this case there is a 
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Log 1( £ 

FIG. 6. The ratio v/u of the four-point couplings u and v as a function of £ for several values 
of s and N = 2. Asymptotically, the ratio converges to v* /u* = 2.9(2) (Ref. [19]). 




FIG. 7. Effective exponents as a function of the rescaled reduced temperature r for several 
values of s and N = 2. 

large downward oscillation. The exponent v& decreases below the asymptotic value and it 
may be even less than 0.5: for s — 1 the minimum value of v e ^ is 0.49 and it is expected 
to further decrease if s increases. These oscillations show how difficult is the determination 
of the critical exponents: extrapolations may provide completely incorrect estimates. The 
effective exponent 7 e ff has a behavior similar to that of u e g. On the other hand, r] e g shows an 
approximately monotonic behavior without detectable oscillations, although the crossover 
effects due to the presence of the 0(4) FP (770(4) = 0.0365(10), Ref. [61]) are clearly visible 
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for s < 1/4. 

B. Crossover behavior in lattice systems 

In Sec. Ill A we computed the FT crossover curves. It is of course of interest to relate 
them to the results obtained in lattice models and in experimental systems. Strictly speak- 
ing, the mapping cannot go beyond the leading correction term appearing in Eq. (3.1) (see 
the discussion in Sec. IV. A of Ref. [52]). In some cases even the leading critical behavior 
cannot be reproduced [54-56]: this happens in the nearest-neighbor Ising model and in the 
lattice self-avoiding walk. However, there are limiting cases in which the FT results exactly 
describe the lattice model: this is the case of the critical crossover limit in weakly coupled 
lattice models and in medium-range models [62-64]. Consider, for instance, a rf-dimensional 
hypercubic lattice and the lattice discretization of the FT Hamiltonian (1.1), 

W= ~^J(x-y)^2<p x , a -<py, a + ^2v(<p x ), (3.6) 

x,y a x 

where the sums over x and y are extended over all lattice points, J(x) is a generic short-range 
coupling, and 

V(<f) =r + ) 2 + § £ • - *«] • (3-7) 

a a ab 

The parameter r is irrelevant and can be made equal to ±1 by changing the normalization 
of the fields. 

The first interesting case corresponds to weakly coupled theories in which r > and 
Uo, Vo 0. Let (3 c (Uo, V ) be the critical point for given U and V and let t be the reduced 
temperature. Then, consider the limit t — > 0, U ,V — > keeping fixed sl = V /U and 
t = t/Ul ,[A ~ d) . In this limit 

X(P,U ,V )U^ fi x F x (als L ), 

t(P,U ,V )U ^^Fz(at,s L ), (3.8) 

where F x (t, s) and F^(t, s) are exactly the FT functions defined in Eq. (3.5). The constants 
/i x , and a can be easily computed by comparing the perturbative expansions (at one loop) 
for the continuum and the lattice model. The additive mass renormalization — it requires 
a nonperturbative matching, see Ref. [63] — also fixes the first terms of the expansion of 
(3 c (Uo, Vo) in powers of Uo and Vq. 

The second interesting case corresponds to medium-range models. In this case we assume 
that the coupling J(x) depends on a parameter p. For instance, one may take 

JJx) = { 1 if M ^ P> (3.9) 

1 otherwise. 

This specific form is not necessary for the discussion that will be presented below, and indeed 
one can consider more general families of couplings, as discussed in Sec. 3 of Ref. [63]. The 
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relevant property is that J p (x) couples all lattice points for p — > oo, i.e., that for p — > oo 
one recovers a mean-field theory. The interaction range is characterized by R defined by 

9 1 V x 2 JJx) 

25 ET^T' (310) 

These models are called medium-range models and admit an interesting scaling limit called 
critical crossover limit [62,64]. If (3 C (R) is the critical temperature as a function of R (here 
[To and Vq are fixed and do not play any role in the limit), then for R — > oo, t = (j3 c (R) — 
P)/P C (R) — > at fixed t = R 2d ^ 4 ~ d h, critical quantities show a scaling behavior. For 
instance, the susceptibility R) and the correlation length scale as 

^^i?)^ 4 - «/«(*)■ (3-11) 

The functions / x (t) and f^(t) are directly related to the crossover functions F x (t,s) and 
F^{t,s) computed in field theory, cf. Eq. (3.5). Indeed, 

fx® = ^x F x( a t,s), 

ft(t)=^Fs(at,s), (3.12) 

where p x , p^, a, and s are nonuniversal constants that depend on the model [63,64]. There- 
fore, the FT crossover functions are expected to describe accurately the crossover behavior 
for large R: in practice, numerical simulations show that for p ps 3 one already obtains a good 
agreement. All constants appearing in Eq. (3.12) can be exactly computed by performing a 
one-loop calculation. The relevant formulae are reported in App. B. 



IV. NUMERICAL RESULTS FROM MONTE CARLO SIMULATIONS 

A. The lattice model 

In order to investigate the existence of the 0(2)®0(2) universality class by numerical 
MC simulations, we considered a simple cubic lattice L 3 and the following Hamiltonian: 

n = -P J2 ■ + i>* ■ Vw) + J2 + €) + (4-i) 

X,fl X 

+a 4 y: m - 1) 2 + - 1) 2 ] + 2A 22 y: 



where <p and ip are two-component real variables. The Hamiltonian Ti describes two identical 
two-component lattice 4 models coupled by an energy-energy term. Note that if A22 = A4 
the symmetry is enlarged to 0(4) and we have the standard four-component lattice 4 model. 
By applying the trasformation 

^1 = ^-, 12 = -— ^, ^ = —jr^ *» = ^/r' (4 - 2) 
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one can easily see that model (4.1) corresponds to the Hamiltonian (3.6) with nearest- 
neighbor coupling J and potential (3.7) with 

U = 12(A 22 + A 4 ), 
Vb = 24(A 22 -A 4 ), 
r = 1 - 2A 4 . (4.3) 

Therefore, model (4.1) is a lattice discretization of the 0(2)®0(2) Hamiltonian (1.1). Ac- 
cording to the FT results presented in Sec. II, continuous transitions in models with A 22 > A 4 
should be controlled by the 0(2)(g>0(2) FP. For A 22 = A 4 the symmetry is enlarged to 0(4) 
and the transition is controlled by the 0(4) FP. If A 22 < A 4 continuous transitions should 
belong to the XY universality class, because the 0(2)<g>0(2) theory has a stable XY FP with 
attraction domain in the region v Q < 0, see, e.g., Ref. [4]. 



B. Monte Carlo simulations 

We present the results of MC simulations for several values of the quartic Hamiltonian 
parameters. We set A 4 = 1 and vary A 22 , considering A 22 = 3, 11/4, 5/2, 9/4, 2, 5/3, and 
7/5. If we define 

they correspond to sl = 1, 14/15, 6/7, 10/13, 2/3, 1/2, and 1/3 respectively. Note that for 
A 22 > At > we have < s L < 2. 

We simulate model (4.1) by using two different types of local moves: (1) a Metropolis 
update in which (p x and ip x are both varied by adding a random term to each component in 
such a way to obtain a 50% acceptance; (2) an 0(4) update [65] in which (p x and ip x are both 
changed keeping fixed the 0(4)-symmetric part of the Hamiltonian, while the 0(4)-breaking 
term 2(A 22 — A 4 ) XLv^V^ is taken into account by performing a standard Metropolis ac- 
ceptance test (the acceptance of this move is rather large, varying from approximately 78% 
for A 22 = 3 to 94% for A 22 = 7/5). In the simulations we use a mixed algorithm in which 
we performed an 0(4) sweep and a standard Metropolis sweep with probability 1/4 and 3/4 
respectively. A rough investigation of the autocorrelation times shows that this is an optimal 
combination. The mixed algorithm is substantially faster (for sl = 1/3 the autocorrelation 
time of the magnetic susceptibility x decreases by approximately a factor of 10) than the 
algorithm in which only the Metropolis update is used. 

We perform a FSS study using lattices with 16 < L < 120 for values of (3 close to f3 c . 
The integrated autocorrelation time r x of the magnetic susceptibility (estimated by using 
the blocking method) increases approximately as r x as cL 2 at (3 C with c 0.2 for A 22 = 7/5 
and c ~ 0.5 for A 22 = 5/3, where the time unit is an update of all spin variables. For larger 
values of A 22 the transition becomes of first order and the dynamics becomes very slow as 
L increases. The large autocorrelation time, i.e. the difficulty of the updating algorithm 
to provide independent configurations, represents the main limitation to the study of the 
critical behavior at f3 c and for large volumes. For each value of (5 we typically performed 
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runs of a few million iterations for the smallest values of L, and of 20-40 million iterations 
for the largest lattice sizes. The total CPU time was approximately 5 CPU years of a single 
64-bit Opteron 246 (2Ghz) processor. 



C. Definitions and notations 

In order to investigate the phase diagram, it is useful to study the FSS of quantities 
related to the energy 

E=yCH) (4.5) 
(V = L 3 is the volume), such as the specific heat 

C es p: «H 2 > - (H) 2 ) , (4.6) 

and the energy cumulant [66] 

S « E, -S' (4 ' 7) 

We also define a quantity M related to the magnetization: 

M=(y/vL*+v%), (4.8) 

where 

/V E yE^' /^ E yE4 ( 4 - 9 ) 

The two-point correlation function G(x) is defined as 

G{x) = (v? • <Px + ipo ■ i/> x )- (4.10) 
The corresponding susceptibility x an d second-moment correlation length £ are given by 

^V GW , ? B ,/ 5(0).- (4 .n) 

V 4 s in 2 (<j mln /2) G (q) 

where G(q) is the Fourier transform of G(x), q = (g m i n ,0,0), and q m \ n = 2n/L. The finite- 
volume definition of £ is not unique. The definition used here has the advantage of a fast 
convergence to the infinite- volume limit [67]. 

In our FSS study we shall consider three RG-invariant ratios [68] 

k = t/L, (4.12) 

H + 



Bi - ' X > (4-13) 



= w-(/v-/v) 2 > (4 14) 

Note that B\ (resp. 5 2 ) is equal to 3/2 (resp. 1/8) at f3 — and to 1 (resp. 0) at f3 — oo. 
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D. First-order transitions: summary of theoretical results 

In the case of a first-order transition the probability distributions of the energy and of 
the magnetization are expected to show a double peak for large values of L. Therefore, as 
a first indication, one usually looks for a double peak in the distribution of the energy and 
of the magnetization. However, as discussed in the literature, see, e.g., Refs. [69,70] and 
references therein, the observation of a double peak in the distribution of the energy for a 
few finite values of L is not sufficient to conclude that the transition is a first-order one. 
For instance, in the two-dimensional Potts model with q = 3 and q = 4 [71,72], double-peak 
distributions are observed for relatively large lattice sizes even if the transition is known 
to be continuous. In order to identify definitely a first-order transition, it is necessary to 
perform a more careful analysis of the large- L scaling properties of the distributions or, 
equivalently, of the specific heat, the energy cumulant, and the Binder cumulants, see, e.g., 
Refs. [66,73,74]. 

The difference of the two maximum values E + and E_ of the energy-density distribution 
gives the latent heat. Alternative estimates of the latent heat can be obtained from the 
lattice-size scaling of the specific heat C and of the energy cumulant Be- According to the 
phenomenological theory [66] of first-order transitions based on the two-Gaussian Ansatz, 
for a lattice of size L there exists a value /9 max of (5 where C has a maximum, C max , and 



Anax-/3c = 0(l/y), 



Cmax — V 



\±1 + 0{1/V) 



(4.15) 



where A h is the (rescaled) latent heat 

A h = E + - E_ 



(4.16) 



Note that, since the temperature parameter f3 is included in the Hamiltonian (4.1), A^ 
should be identified with the dimensionless ratio between the latent heat and the critical 
temperature. We recall that in the case of a continuous transition one expects 



Anax — Pc ~ a L i C max ~ bL a ^ U + C. 



(4.17) 



The energy cumulant Be can also be used to identify first-order transitions. Indeed, a careful 
analysis [66] shows that there is a value P min where B E has a minimum, B E ,mm, and which 
is related to the latent heat. The phenomenological theory gives [75] 



Pmin-Pc = 0(l/V), 



B 



i?,min 



2 



+o(i/n 



(4.18) 



where 



E+ — E— 
VE^EZ' 



(4.19) 



In continuous transitions A = — there is only one peak in the energy distribution — and the 
infinite-volume limit of B E ,m.m is trivial: Hindoo Be, mm — 2/3. 
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As discussed in Ref. [73], the distribution of the order parameter is also expected to show 
two peaks at M + and M_, M_ < M + , with M_ — > as L — > oo since in the high-temperature 
phase there is no spontaneous magnetization. The phenomenological theory predicts that 
the Binder parameter can still be used to identify the critical point [the analysis shows 
that /3 C ross (Li, 1/2 ) — Pc ~ min(Li, L 2 )~ 2d , where j3 CTOSS (Li, L 2 ) is the value of /3 at which 
B 1 (L 1 ) = Bi(L 2 )]. Moreover, it predicts 

dBjdpy^ » cL d , (4.20) 

for sufficiently large lattice sizes, where d is the space dimension. More generally, close to 
P c the phenomenological theory predicts [73,74] 

Bi(P, L) ~ - Pc)L d \. (4.21) 

Such a relation is valid only sufficiently close to P c since the scaling function f(x) diverges for 
x = a: pea k where x peSb k is related to the position of the peak present in B x for L — > oo. Ref. [73] 
indeed shows that Bi(P, L) at fixed L has a maximum £? lmax (L) at P — P pea k(L) < P c with 

B hmax ~ L d , /3 peak - & « c p L- d . (4.22) 

Thus, for (P — P c )L d close to (/3 pca k — /3 c )-^ d = a^pcak = c P the scaling behavior (4.21) breaks 
down [/(xpeak) diverges] and subleading terms of order L~ d must be included. In the region 
in which Eq. (4.21) holds, we have dB x jdp « L d /'[(/5 - /9 c )L d ]. By using Eq. (4.21) we can 
express (P — P c )L d in terms of B\ obtaining the relation 

^ = Ld fW, (4.23) 

with a suitable scaling function f(x). In practice this means that L~ d dB\j dp converges to 
a universal function of B\ as L — > oo. 

It should be noted that the phenomenological theory has been developed for systems with 
discrete symmetry group, for instance for the Potts model [66,73]. In our case the symmetry 
is continuous; thus, one may wonder if it can be also applied to the present case. The 
numerical results that we will present below in Sec. IV E show that no changes are needed 
and that all predictions hold irrespective of the symmetry group. This can be understood on 
the basis of a simple argument. Imagine we introduce a magnetic field H in the model. The 
first-order transition should be robust with respect to H, since we expect here the transition 
to be temperature driven. In other words, we expect the behavior to be unchanged if H is 
switched on, as long as H is small. For H ^ we have a discrete system, thus all previous 
scalings apply. If the behavior is continuous in H, all results also apply for H = 0. Note 
that, for H ^ 0, B 2 is expected to be noncritical, since B 2 vanishes in a system magnetized 
in a specific fixed direction. Thus, B 2 should have no discontinuity at the transition and 
its derivative at P — P c should be finite, at variance with the behavior of B±, cf. Eq. (4.20). 
This fact will also be verified by the MC results that we shall present in Sec. IV E. 
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FIG. 8. Hystograms of the energy density (above) and of the magnetization (below) for A22 = 2 
and lattice sizes L = 80, 100, 120. For each L we report the data for the value of (3 at which the 
two peaks have approximately the same height. These results have been obtained by reweighting 
the MC data at /3 = 0.89170, 0.89175, 0.81980. 

E. The phase diagram for A22 > A4 = I 

In this subsection we investigate the phase diagram of the lattice model (4.1) for A22 > 
A4 = 1 with the purpose of identifying the regions in which the model shows a first-order or 
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FIG. 9. Specific heat for A22 = 5/2 and several values of L. 

a second-order phase transition. 

The distributions of E and M show two peaks for A22 > 2 when L is large enough. For 
A 2 2 = 3 two peaks are already observed for L = 16, while for A 2 2 = 2 two peaks are observed 
only for L > 100, cf. Fig. 8. Thus, the model with A4 = 1 has apparently a relatively strong 
first-order transition for A22 > 3 that gradually weakens as A22 is decreased. In order to 
check that we are really in the presence of a first-order transition, we check the scaling 
behavior of the different observables for L — ► 00. The predictions for C and Be are well 
verified. For instance, in Fig. 9 we show C/V for A 2 2 = 5/2 for several values of L. In 
agreement with Eq. (4.15) C max /V has a finite limit as L — ► 00. In Fig. 10 we show the 
energy cumulant Be for the same value of A 2 2', Be, m in is different from 2/3, confirming the 
first-order nature of the transition. From Be, mm we can determine the latent heat. For each 
L we determine A^, which is obtained from BE,mm(L) by using Eq. (4.18) and neglecting the 
1/V corrections. The latent heat is obtained by extrapolating assuming 1/V corrections. 

Estimates of Al for several values of A22 and L are reported in Fig. 11. They show the 
expected 1/V behavior and allow a precise determination of the latent heat in the infinite- 
volume limit. The results of the extrapolations are reported in Table II. They are in perfect 
agreement with those obtained from the maximum of the specific heat, using Eq. (4.15), and 
from the position of the peaks in the energy distributions. 

We have also performed MC simulations for A22 = 5/3 and 7/5 on lattices of size L < 120, 
without observing evidences for first-order transitions. The hystograms of E and M do 
not show any evidence of two peaks and are not significantly broad, as for instance the 
distribution of E for L = 80 and A 2 2 = 2, cf. Fig. 8, which would indicate the onset of 
two peaks. In the cases A22 = 5/3, 7/5, C and Be on lattices of size L < 120 do not show 
the behavior predicted by the phenomenological theory of first-order transitions [66]. In 
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FIG. 10. Energy cumulant Be for A22 = 5/2 and several values of L. 



particular, B Erain continuously increases towards 2/3 and we can only put upper bounds on 
A h . For A 2 2 = 5/3 we obtain for instance < 0.005. A more stringent bound is suggested 
by the width of the energy distribution at /3 C , i.e. A^ < 0.003. 

Let us now discuss the behavior of the variables related to the magnetization. In par- 
ticular, we focus on the derivatives with respect to (3 of B\ and B 2 . Predictions (4.22) and 
(4.23) are well verified by our data with A 2 2 > 2. We observe the presence of a peak in 
B\ at fixed L that becomes sharper as L increases and we also verify that far from this 
peak Eq. (4.23) holds. Of course, corrections increase as A22 decreases, as expected. On the 
other hand, B\ is monotonically decreasing with /3 at fixed L for A22 = 7/5 and A22 = 5/3, 
providing no evidence that the transition is of first order for these two values of A 2 2- 

We have repeated the same analysis for B 2 . Its behavior looks quite different with respect 



TABLE II. For the model with A4 = 1 and several values of A22 we report sl, /3c; the average 
energy density E c = (E + + EJ)/2, A, cf. Eq. (4.19), and the latent heat A/j. For (3 C and E c all 
reported digits are exact. 



A22 sl (3c Ec A A h 

3 1 0.8733 1.98 0.0607(6) 0.1198(12) 

11/4 14/15 0.8627 1.98 0.0396(12) 0.0783(24) 

5/2 6/7 0.8504 1.97 0.0239(4) 0.0471(8) 

9/4 10/13 0.8364 1.96 0.0120(4) 0.0235(8) 

2 2/3 0.8198 1.94 0.0048(2) 0.0093(4) 

5/3 1/2 0.7927 1.92 < 0.0015 < 0.003 
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FIG. 11. Values of for A22 = 2, 9/4, 5/2, 11/4, 3, as obtained from BE, m m by using 
Eq. (4.18). The lines represent the infinite- volume extrapolations assuming a 1/V asymptotic 
correction. 

to B\. First of all, we observe a monotonic behavior without peaks for all the considered 
values of A22, including the largest ones for which the peak in B\ is rather sharp. Moreover, 
the data are reasonably well described by assuming 



suggesting that B 2 does not have a jump at the transition, as expected on the basis of the 
argument presented at the end of Sec. IV D. A similar analysis can also be performed for 
1% for which we have no prediction. Our data are roughly consistent with a behavior of the 
form 



with 3.2 < 9 < 3.7. It is quite difficult to intepret such a value of 9 that could well be d or 
d + 1 depending on the size of the corrections. 

The results for A/, reported in Fig. 12 suggest therefore a phase diagram characterized 
by a line of first-order transitions extending from large values of A22 down to A 22 where A^ 
vanishes. In order to compute A 22 , we should first discuss the expected behavior of close 
to the tricritical point A 22 (see, e.g., Ref. [76]). We consider a generic model depending on 
two parameters f3 and g with a tricritical point at g*, (3* = (3 c (g*). The critical behavior 
can be parametrized in terms of two linear scaling fields u\ and U2 with RG dimensions y\ 



dB 2 



f(Bx), 



(4.24) 



d/3 



dp 



L 9 f(Bi) 



(4.25) 
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FIG. 12. Latent heat A^ versus ^22- The line corresponds to the best fit, 
A^j = 0.049(^22 — 1.52) 2 ' 29 . We also report the estimate of A 22 with the corresponding error: 
A* 22 = 1.52(6). 

and y 2 , satisfying y\ > y 2 . In the absence of any symmetry, the linear scaling fields are 
combinations of (g — g*) and of (/? — /?*). The scaling field u\ is completely defined, while 
u 2 can be arbitrarily chosen as long as it is independent of U\. Therefore, we write 

Ui = axifl - (3*) + a 2 (g - g*), 

u 2 = g-g*. (4.26) 

The first-order transition line is characterized by the equation u 2 = ciuf 2 1 with, by a 
proper choice of the scaling field, U\ > 0. Analogously, the second-order transition line is 
given by u 2 = c 2 (sui) y2 ^ yi where s may be either 1 or —1 and |ci| 7^ | C2 1 • Note that both 
lines are tangent to U\ = 0, but that the nonanalytic deviations are parametrized differently. 
The free energy can be written as 

F w F reg ((3,g) + \u 2 \ d ^F sing>± ( Ul \u 2 \-y^), (4.27) 

where i'kng^ (resp. F sing _) applies to the case u 2 > (resp. u 2 < 0), and F reg (/3,g) is 
a regular function. By hypothesis, assuming c\ > without loss of generality, F singi+ (x) 
is continuous but has a discontinuous derivative at the first-order transition line, i.e. for 
x = c - yi/y2 . It follows 

E + -E^^ a 2 \u 2 \ {d -^ )/y2 /\F' s ^ + ~\g- g*]^^. (4.28) 

In the model we consider g = A 22 , so that we predict that sufficiently close to the tricritical 
point A 22 
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A, = A (A 22 - A* 22 f, 



(4.29) 



where 9 is an exponent defined by the tricritical theory. If we fit our data for the latent heat 
with this expression we obtain 

A* 22 = 1.52(6), A = 0.049(7), 9 = 2.29(15), (4.30) 

with a x 2 P er degree of freedom (DOF) of 0.24. The fit is stable with respect to the 
number of points that are included: If the point with A22 = 3 is not included we obtain 
A 22 = 1.50(12) and 9 = 2.35(36), in full agreement with the result reported above. We have 
also investigated the possibility that A 22 = A4 = 1, which would imply the absence of critical 
transitions in the 0(2)(g>0(2) universality class. A fit of the data for the latent heat of the 
form A h = A (A 22 -l) e gives 9 = 3.49(4) with xV D0F ~ 9 - If the P oint with ^22 = 3 is not 
included we obtain 9 = 3.87(4) (x 2 /DOF ~ 11), while if the two points with largest A22 are 
discarded we obtain 9 = 4.00(3) (% 2 /DOF « 9). These fits are significantly worse than that 
presented above. Moreover, the estimate of 9 does not agree with the theoretical prediction 
that can be obtained assuming the tricritical point to be at A 2 2 = 1- Indeed, if A 22 = 1 the 
tricritical theory coincides with the 0(4) theory. There are two relevant perturbations, the 
thermal perturbation with RG dimension 1/v and the perturbation that breaks the 0(4) 
symmetry down to 0(2)£g>0(2), which is associated in field theory with the spin-4 quartic 
operator, with RG dimension <p/v. For the 0(4) universality class, a = —0.247(6) (Ref. [61]) 
and = 0.08(1) (Ref. [48]). Since < 1 in this case, we have yi = 1/u and y 2 = so 
that we should have [77] 

A, ~ (A22 - l) {l - a)/4> ~ (A22 - 1) 16±2 , (4.31) 

which is not compatible with the MC results. Thus, our results are incompatible with 
A 22 = 1, unless the limiting behavior sets in for even smaller values of A^. This seems 
unlikely, since for A22 < 3 we are already in a region in which A^ <C 1. Thus, our simulations 
provide a nice evidence in favor of the existence of a tricritical point at A 22 > 1 separating 
the first-order transition line from a continuous transition line for 1 < A 2 2 < A 22 . Notice 
that a phase transition is always expected due to the qualitative difference of the minimum 
of the free energy in the high- and low-temperature regions. 

The estimate (4.30) implies that for A22 = 5/3 one should observe a first-order transition. 
As we have already discussed, our data do not show any evidence for that. This is hardly 
surprising, since the fit of the estimates of A h predicts A h = 0.8(6) • 10~ 3 for A 2 2 = 5/3, 
which is smaller than the bound A h < 3 ■ 10~ 3 obtained on lattices L < 100. A posteriori, 
however, one can convince oneself of the first-order nature of the transition by considering 
the estimated values of v. For instance, the derivative of the Binder parameter B\ with 
respect to (3 is expected to scale at the critical point as B[ ~ L d for a first-order transition. 
In our analyses of the data with A 2 2 > 2, we find that, if one estimates B[ with B[ ~ L e , 
the effective 9 rapidly increases towards d as the considered set of sizes increases, even when 
a double peak in the distributions is not yet evident. For A22 = 5/3 if we use the data 
with L > 60 we obtain 1/9 = 0.47(2), while including only the data with L > 80 gives 
1/9 = 0.42(4): the exponent 9 is increasing towards the expected value 9 — d — 3. At a 
second-order transition one expects B[ ~ L 1 ^, and thus the previous results would give a 
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FIG. 13. The inverse critical temperature j3 c versus A22- The vertical lines indicate the estimate 
of A 22 with its error. The data for A22 > A 22 = 1.52(6) are indicated by circles, while the ones for 
A22 < A 22 are indicated by a square. Errors are much smaller than the size of the symbols. The 
dashed line shows the linear fit of the three data for A22 > A 22 that are closest to A 22 . 



quite small value for z/, definitely incompatible with the FT predictions. Therefore, even if 
we do not have a direct evidence, the data for A22 = 5/3 are better explained by a very 
weak first-order transition than by a second-order one. 

The presence of a tricritical point at A 22 = 1.52(6) is further supported by the results for 
j3 c reported in Table II and plotted in Fig. 13 versus A 2 2 — A*. Equation (4.26) implies that, 
sufficiently close to the tricritical point A 22 , the values of f3 c along the first-order transition 
line behave as 



-^{A 2 2 



a;,) + - 



a 1 



-(A 22 

Cl 



1 yi/y2 



A*22) 



(4.32) 



while, along the second-order transition line, we have similarly 



f3 c -(3:^--(A 2 2-A* 22 ) + — 
ai sa± 



c 2 



■{A 22 -A* 22 ) 



1 yi/yz 



(4.33) 



Equations (4.32) and (4.33) strictly hold only if yi/y2 < 2, otherwise one should also include 
additional integer powers (A22 — A 22 ) n , n < yi/i/2, with coefficients that are identical for the 
first-order and second-order transition line. 

The linear dependence of (5 C on A22 for A 2 2 > A 22 is clearly observed in Fig. 13. A linear 
fit of the three points that are closest to A 22 and satisfy A 2 2 > A 22 gives a 2 /ai ps —0.068 
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with a reasonable \ 2 ■ Note that the deviations from a straight line are very small, indicating 
that the corrections, and in particular, the nonanalytic ones, are tiny. Fig. 13 also shows the 
value of f3 c for A 22 = 7/5 < A 22 , i.e. f3 c = 0.76615(15) which will be determined in the next 
subsection, and for the 0(4) 4 model obtained by setting A 22 — A = 1, which is given by 
P c = 0.7154(1) [78]. Both values, and in particular the one for A 22 = 7/5, differ significantly 
from the linear extrapolation of the data for A 22 > A 22 . This behavior of f3 c as a function 
of A 22 is naturally explained by the presence of a tricritical point at A 22 = 1.52(6), and it is 
another evidence of the fact that the point A 22 = 7/5 belongs to the second-order transition 
region. Indeed, in this scenario the nonanalytic corrections are different (|c 2 | ^ |ci|) from 
those observed in the first-order transition region. Thus, the observed value for A 22 = 7/5 
can be explained by the presence of nonanalytic corrections to the linear behavior that are 
substantially larger (\c 2 \ <C |ci|) than those observed on the side of the first-order transition 
line. The data shown in Fig. 13 can hardly be explained by assuming A 22 = 1, i.e. a first- 
order transition line extending down to the 0(4) point A 22 = A 4 = 1. Indeed, in this case 
the linear behavior should extend down to the 0(4) point, but this is clearly contradicted by 
the values of f3 c obtained for A 22 = 7/5 and A 22 = 1. Note also that the same RG arguments 
leading to Eqs. (4.26), (4.27) and (4.32) tell us that the values of j3 c along the continuous 
transition line must approach linearly the 0(4) point A 22 = A, but with a slope different 
from the one observed at the tricritical point A\ 2 = 1.56(2). 

In conclusion we have shown that, in the quartic parameter space with A 22 > At, there 
is a region in which the transition is continuous and therefore belongs to the 0(2)00(2) 
universality class controlled by the 0(2)<g>0(2) FP found in the FT approach of Sec. II. 

Finally, we should note that in the FT studies the basin of attraction of the 0(2)<g>0(2) 
FP includes all theories with s < s* and s* > 1, while here s* L = 0.41(4). As we have already 
discussed before, in spite of the similar definition, sl should not be identified with s and we 
have sl — > s only for a weakly coupled theory, i.e. for A — > and A 22 — > 0. Here, A = 1, 
so that we are far from this limiting case. In any case, the FT results imply that s* L increases 
if A decreases. The critical value s* L is also expected to increase if longer-range interactions 
are added. Indeed, as shown in App. B, we expect 0.71 < s* L < 1.23 for medium-range 
models with A^ — l. 



In this Section we analyze the results for A 22 = 7/5. On the basis of the analysis 
presented in the previous Section, in this case the transition should be of second order. In 
the FSS limit we expect the following behavior 



where R is any RG invariant quantity (we will take R to be Bi, B 2 , or l^). These scaling 
forms are valid for (3 — > (3 C) L — > 00 at fixed argument (f3 — f3 c )L l / y . From Eq. (4.34) we 
obtain moreover 



F. Determination of the critical quantities for A 22 = 7/5 



R^MW-^L 1 ^}, 



(4.34) 
(4.35) 
(4.36) 
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TABLE III. Analysis of the RG-invariant quantities R, R = B\, B 2 , and l^, by using Eq. (4.34). 
We report u, C , and R* = R{f3 c ). 



7-) 

R 


T 

^min 


X /DOF 


V 


Pc 


R 


Bi 


20 


78/70 


0.728(4) 


0.766284(8) 


1.1111(6) 




40 


63/54 


0.716(7) 


0.766278(11) 


1.1121(9) 




60 


40/38 


0.711(14) 


0.766281(17) 


1.1121(16) 




80 


24/20 


0.665(27) 


0.766260(25) 


1.1148(30) 


B 2 


20 


305/69 


0.662(5) 


0.765890(7) 


0.0643(2) 






Q9 KA 


u.uoy yv j 


U. t Uoao^ylV ) 






60 


51/38 


0.650(17) 


0.766044(15) 


0.0580(7) 




80 


25/21 


0.587(33) 


0.766038(19) 


0.0585(9) 


k 


20 


117/71 


0.687(2) 


0.766153(3) 


0.5692(6) 




40 


61/55 


0.692(2) 


0.766168(5) 


0.5726(10) 




60 


29/37 


0.694(5) 


0.766181(8) 


0.5770(19) 




80 


9/18 


0.673(13) 


0.766168(12) 


0.5741(39) 



^^L^f' R [(P-p c )L^]. (4.37) 

Eqs. (4.35), (4.36), and (4.37) can be used to determine 77, f3/v, and v. It is also possible 
to avoid the use of the two unknown quantities, f3 c and v. We use Eq. (4.34) to express 
(P — fi c )L x l v in terms of R, and rewrite all equations as 

X~L 2 -*g x , R (R), (4.38) 

M*L-V»g M , R (R), (4.39) 
dR 1 



dp 



L^g RuR2 (R 2 ), (4.40) 



where R is a RG-invariant quantity. 

Let us begin by performing a direct analysis of B\, B2, and l^. We fit our results 
corresponding to 20 < L < 120 and 0.755 < P < 0.767 (see Fig. 14) by using Eq. (4.34). 
For this purpose we must somehow parametrize the scaling function f R (x). We use a simple 
polynomial expression, writing 

v 

f R ( x) = Y^a n x n . (4.41) 

n=0 

The order p is chosen in the following way. For a given set of data we perform a nonlinear 
fit, increasing each time p until the x 2 changes by approximately 1 by going from p to 
p+1. Of course, one should also worry about scaling corrections and crossover effects. In 
order to detect them we perform the fit several times, each time including only the data 
corresponding to lattice sizes L larger than some value L min . The corresponding results are 
reported in Table III. 

The quality of the fits of B\ and is reasonably good (x 2 /DOF 1), while for B 2 one 
should certainly discard the results with L min = 20 and 40. As far as the estimates of v are 
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FIG. 14. Plots of h, B\, and B2 for A22 = 7/5. The lines are drawn to guide the eye. 
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concerned, we observe in all cases a systematic drift. The analyses of B\ and give first 
v ~ 0.69 - 0.71 and then, by using only the data with L min = 80, one obtains v ~ 0.66-0.67. 
On the other hand, fits of B 2 give first v ps 0.66 and then, for L min = 80, v ~ 0.59, although 
with a large error of ±0.03. Clearly the data are affected by large scaling corrections and, 
apparently, even on these relatively large lattices one is not able to obtain a precise estimate 
of u, but only an upper bound v < 0.67. In any case, if we assume that the observed 
discrepancies give a reasonable estimate of the systematic error, the results with L min = 80 
give 



which includes the estimates from Bi, B 2 , and with their errors. This is fully compatible 
with the FT estimates. 

The results obtained from the analysis can be interpreted in terms of a crossover due to 
the presence of a nearby 0(4) FP. Indeed, the observed behavior rensembles quite closely 
what is observed in field theory for, say, s = 1/4. The effective exponent v is first close to the 
0(4) value and then decreases towards its asymptotic value. This interpretation is somehow 
supported by the observed values of B* x and Z|, where R* = R(f3 c ). They are close to the 
corresponding 0(4) ones [61]: /| « 0.547, B{ = 1.0945(2). Only for B 2 do we observe a 
relatively large difference since in the 0(4) case B* = ±B{ = 0.09121(2). This may explain 
why the estimates of v from B 2 are those that most differ from the 0(4) ones. Of course, 
one may not exclude that the asymptotic values B{ and /| are close to the 0(4) estimates. 

As a check we have also considered the derivatives of B 1 , B 2 , and 1% and we have used 
Eq. (4.40). We find that in all cases the best fit (smallest x 2 /DOF) is obtained by taking 
R 1 = R 2 with results that are (not surprisingly) fully compatible with those obtained in 
the previous analysis. Fits with R\ ^ R 2 are somewhat worse but always show the same 
pattern. While for small values of L m [ n , v varies between 0.64 and 0.70, depending on the 
choice of R\ and R 2 , for L min = 80 the analyses indicate a smaller value, fully compatible 
with the result reported above. 

The analyses of B 1 , B 2 , and 1% also provide estimates of f3 c . There is a clear upward trend 
in the estimates obtained from B 2 , as it can be also understood from Fig. 14: the value at 
which the lines B 2 (f3, L) at fixed L cross moves significantly towards higher values of (3 as L 
increases. The estimates obtained from B\ and 1% are apparently stable, but not compatible 
within the tiny statistical errors indicating that there are strong (compared to the statistical 
errors) crossover effects, as observed for v. If we assume that the discrepancies among the 
estimates obtained from the analyses of the three RG-invariant quantities give a reasonable 
estimate of the systematic error, an estimate of f3 c that includes all results is 



We finally compute 77 and (3/v from the analysis of x an d M. We have performed the 
analyses by using Eqs. (4.38) and (4.39). The analyses using R — 1% are well behaved and 
show little dependence on L min , leading to the estimates 



v = 0.63(7), 



(4.42) 



P c = 0.76615(15). 



(4.43) 



r] = 0.045(10), 



0.0525(10). 



(4.44) 
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The corresponding scaling plots are reported in Fig. 15. On the other hand we observe 
systematic deviations if we use R — B\ or B 2 . The goodness of the fit, x 2 /DOF, is a factor 
of ten larger than for the analysis with R = 1%, and the estimates vary significantly, between 
-0.1 and 0.2 (77) and 0.45 and 0.6 (/3/V)- 

The results (4.44) satisfy the scaling relation 2/3/u = 1 + rj quite precisely They can be 
compared with the FT results: rj = 0.09(1) (MZM scheme, Ref. [19]) and 77 = 0.09(4) (MS 
scheme). Although larger, the MS estimate is compatible within error bars. A discrepancy 
is observed for the MZM result, whose error might have been underestimated (after all, the 
MZM series are not Borel summable). Of course, it is also possible that scaling corrections 
play an important role, as it is the case in the analyses of Bi, B 2 , and 1%. In this case, 
we expect the MC result to be influenced by the presence of the nearby 0(4) FP. Such 
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an interpretation is supported by the fact that the estimated rj is close to the 0(4) result 
rj = 0.0365(10) (Ref. [61]). 

Finally, we analyzed the specific heat by using 

C = L a ' u f(R), (4.45) 

which is valid as long as a > 0. We do not expect this fit to be very precise since we are 
neglecting the analytic contribution that gives rise to corrections of order L~ a / U , which are 
expected to be sizeable since a is small. All fits have a poor % 2 . Only the fit with R = 1$ 
and L m i„ = 80 has x 2 /DOF of order one. Fits using 1% give estimates of a/u that decrease 
as L m i„ increases, varying between 0.20 and 0.16. Fits using B\ show the same decreasing 
trend with 0.19 < a/u < 0.20. Fits using B 2 — they have a very large x 2 , X 2 /DOF = 14 for 
Lmin = 80 — show a more erratic behavior with L min and give 0.14 < a/u < 0.18. We quote 
as final result that obtained by using and L min = 80: 

^ = 0.16(3). (4.46) 

The error has been chosen such that it includes all estimates with L min = 80. Using hyper- 
scaling Eq. (4.46) gives u = 0.63(1), in good agreement with the result reported above, and 
a = 0.10(2). A scaling plot is reported in Fig. 16. 

The very large x 2 /DOF of the above-reported analysis is probably due the fact that the 
analytic contribution to the specific heat is neglected. We have thus performed a second set 
of analyses using 

C = L a ^f 1 (R) + f 2 (f3) (4.47) 

and R — £/L. As before, we have used polynomials for fi(x) and f2(x). Fifth-order 
polynomials allow us to obtain x 2 /DOF ps 1 even for L min = 20. We obtain a/u = —0.10(6), 
0.01(9), and 0.16(25) for L min = 20, 40, 60 respectively. The results for the largest L min are 
compatible with the estimates reported above. However, the very large errors indicate that 
our data are not precise enough and that our set of values of L is too small to disentangle 
the analytic background from the singular behavior. This probably means that the error 
reported in Eq. (4.46) should not be taken too seriously and is most likely underestimated. 



V. CONCLUSIONS 

In this paper we investigated the critical behavior of three-dimensional models with 
symmetry 0(2)<g>0(iV) described by the FT Hamiltonian (1.1) in the case i>o > 0, which 
corresponds to the symmetry-breaking pattern 0(2)®0(iV)— >0(2)®0(iV — 2) (the case vq < 
has been discussed in detail in Ref. [11]). 

First, we considered the FT perturbative approach. The analysis of the five-loop series 
in the MS scheme without e expansion provides strong evidence for the presence of a stable 
FP with v > for N = 2 and N = 3, and therefore for the existence of the corresponding 
three-dimensional 0(2)(g>0(iV) universality classes. This result confirms the conclusions of 
Ref. [19], in which a stable FP with v > was found in the three-dimensional MZM scheme. 
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Note that these FT perturbative analyses disagree with the conclusions of Refs. [30,31,5], in 
which no FP was found by using a nonperturbative RG approach. Moreover, the MS scheme 
without e expansion allowed us to obtain fixed-dimension results at any d. We recovered the 
e-expansion results sufficiently close to four dimensions and obtained a full picture of the 
fate of the different FPs as d varies from four to three dimensions. 

In order to confirm the existence of a new universality class we have performed a MC 
study of a lattice discretization of Hamiltonian (1.1) for N = 2. The purpose is to identify 
a parameter region in which the transition is of second order with the expect symmetry- 
breaking pattern, 0(2)®0(2)^0(2). A detailed analysis of the critical behavior of the model 
(1.2) for A4 = 1 and A22 > A4 shows the following phase diagram. For A22 > A\ 2 = 1.52(6) 
the model has a first-order transition. Such a transition is easily identified for A22 > 3 — the 
energy and the magnetization show two peaks already for L « 10 with reduced latent heat 
AE/T C larger than 0.1. The first-order transition becomes weaker as A 2 2 decreases, the latent 
heat vanishing at the tricritical point A 2 2 = A* 22 . For 1 = v4 4 < A 2 2 < A* 22 the transition is 
continuous. Its critical behavior is expected to belong to the three-dimensional 0(2)(g>0(2) 
universality class and therefore to be controlled by the stable FP of the 0(2)® 0(2) theory 
found within the FT methods. We have also performed simulations of the lattice model for 
^22 = 7/5 in order to identify the critical behavior. All results are definitely compatible with 
the expected behavior at a second-order transition. We are however unable to provide precise 
estimates of the critical exponents since we observe strong crossover effects, probably due to 
the presence of the nearby 0(4) FP. The effective exponents computed in the MC simulation 
rensemble those observed in the FT model for small s = Vq/uq, see Fig. 7, indicating that 
crossover effects play an important role in these systems and make difficult, both numerically 
and experimentally, a precise determination of the asymptotic critical behavior. 

The FT Hamiltonian (1.1) is supposed to describe the critical behavior of STA's and of 
helimagnets [3]. Inelastic neutron-scattering experiments show that STA's can be modeled 
by three-component spin variables s associated with each site of a stacked triangular lattice 
and by the Hamiltonian 
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(vw) xy {vw) z v 



(5.1) 



The first sum is over nearest-neighbor pairs within the triangular layers (xy planes) with 
an antiferromagnetic coupling J\\ < 0, the second one is over orthogonal interlayer nearest 
neighbors. If the uniaxial term is positive, one has an effective two-component theory. 
Numerically, Hamiltonian (5.1) has been much studied in the limiting cases D = +oo and 
D = 0. In the first case the spins are confined to a plane, i.e., one is effectively considering 
XY spins. There is now evidence that XY STA's undergo first-order transitions, at least 
for | J±/J\\\ not too small. Indeed, first-order behavior has been observed for J±/ J\\ = —3/4 
(Ref. [33]), JJJ\\ = -1 (Ref. [34]), and J ± /J\\ = -10 (Ref. [32]). The numerical results 
indicate that the first-order transition becomes stronger as \J±/J\\\ increases, and thus we 
can conclude that all these models have first-order transitions, at least for |Jj_/J||| > 3/4. 
It should however be noted that the latent heat is very small. For J±/ J\\ = —3/4 and 
J±/J\\ = -1, AE/T C ~ 7 • 1(T 3 [33,34]. This means that small modifications of the lattice 
model may turn the first-order transition into a second-order one. In particular, it is not 
clear whether, on the basis of these numerical simulations of the XY STA's, we should 
expect first-order transitions also for experimental easy-plane systems, which do not satisfy 
the condition D 3> J\\, \ J±\. For instance, in the case of CsMnBr 3 we have [1] J\\ « 0.0018 
meV, Jj_ = 0.88 meV, and D = 0.013 meV, while in other compounds like CsVX 3 (X = 
CI, Br, I) one observes [1] D « Jy <C J±. 0(2)®0(2) critical behavior is also expected 
in easy-axis materials (D < 0) in the presence of a (large) magnetic field along the easy 
direction. 

Experiments on STA's favor a second-order transition, although the estimates of (3 do 
not satisfy the inequality (3 > v/2 that is expected on the basis of unitarity [31,5]. Of course, 
this could be explained by the presence of a weak first-order transition. A second possibility 
is that the experimental systems are in the basin of attraction of the stable FP, but close 
to the boundary of the stability region. If this is the case, on the basis of the FT crossover 
analysis, one expects strong crossover effects; see, for example, the results presented in 
Fig. 7 for s ~ 1. Thus, the exponents (3 and v that are measured experimentally may well 
differ from their asymptotic value. For what concerns helimagnets, their critical behavior 
is somewhat different from that observed in STA's. The exponent (3 is always very close to 
the 0(4) value, while a varies between 0.1 and 0.3. These results strongly remind our MC 
ones with A 2 2 = 7/5. In that case, (3/v was close to the 0(4) value and a w 0.10. Thus, the 
helimagnetic results can be explained by the presence of the nearby 0(4) FP that controls 
the critical behavior for \t\ > 10~ 3 . 

Finally, we would like to conclude with some remarks on the experimental relevance of 
the 0(2)®0(3) universality class. Such a critical behavior is expected in some easy-axis ma- 
terials that have a small uniaxial anisotropy [1], for instance in RbNiCl 3 , VC1 2 , and VBr 2 . 
However, the reduced-temperature region in which 0(2)®0(3) behavior might be observed 
is usually very small, i.e. for t > 10~ 2 , because these systems are expected to crossover 
to an XY critical behavior for t < 10~ 2 [1,3]. Therefore, the asymptotic 0(2)®0(3) criti- 
cal behavior can be hardly observed in these materials and significant differences between 
theoretical predictions and experimental results should not be unexpected. As argued in 
Ref. [79], and usually assumed in the literature, the 0(2)®0(3) critical behavior should also 
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be experimentally realized in easy-axis STA's, such as CSMCI3 and CsNiBr3, at the multi- 
critical point observed in the presence of an external magnetic field along the easy axis, or 
at the critical concentration of mixtures of easy-axis and easy-plane materials, for instance 
in CsMn(Br :E I 1 _ :E )3 [80]. We note that the identification of the multicritical point with the 
0(2) (8)0 (3) universality class is not obvious and should be theoretically analyzed. The criti- 
cal behavior at the multicritical point in a magnetic field should be described by the stable FP 
of the most general LGW theory with symmetry 0(2)<g)[Z 2 ©0(2)] [81]. In Ref. [79] only the 
quadratic terms have been considered and discussed, but the relevant LGW Hamiltonian has 
also additional quartic terms beside those appearing in the 0(2)® 0(3) Hamiltonian. As a 
consequence, the 0(2)®0(3) FP, describing a critical behavior with an enlarged 0(2)®0(3) 
symmetry, determines the asymptotic behavior at the multicritical point only if it remains 
stable with respect to the additional quartic terms breaking 0(2)(g>0(3) to 0(2)(g>[Z 2 ©0(2)]. 
The critical behavior at the multicritical point is determined by the stable FP of the RG 
flow of the complete LGW theory with symmetry 0(2)®[Z 2 ©0(2)]. This issue was recently 
investigated in Ref. [81] by a FT analysis based on five-loop calculations within the MS and 
MZM schemes. Unfortunately, this study was unable to establish the stability properties of 
the 0(2) ©0(3) FP. In any case, it did not provide evidence for any other stable FP. Thus, 
on the basis of these FT results, the transition at the multicritical point is expected to be 
either continuous and controlled by the 0(2) ©0(3) fixed point or to be of first order. Similar 
arguments can be applied to the multicritical point in mixtures of easy-axis and easy-plane 
materials, such as CsMn(Br x Ii_ a; )3. We believe that the identification of the multicritical 
behavior with the 0(2)©0(3) universality class is even more questionable in this case, since, 
beside the additional quartic terms considered above, there are other perturbations related 
to the quenched randomness. 
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APPENDIX A: THE FIVE-LOOP SERIES OF THE MS SCHEME 

In the MS scheme [37] one sets 

0= [Z^v)] 1 ' 2 ^ (Al) 
u = A d fi € Z u (u,v), 
v = A d fi e Z v (u,v), 

where the renormalization functions Z^, Z u , and Z v are determined from the divergent 
part of the two- and four-point one-particle irreducible correlation functions computed in 
dimensional regularization. They are normalized so that Z^u^v) ~ 1, Z u (u,v) ~ u, and 
Z v (u,v) ~ v at tree level. Here Ad is a (i-dependent constant given by Ad = 2 d ~ 1 ir d / 2 T(d/2). 
Moreover, one defines a mass renormalization constant Z t (f,g) by requiring ZtT^ 1 ' 2 ^ to be 
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finite when expressed in terms of u and v. Here I^ 1 ' 2 ) is the one-particle irreducible two-point 



function with an insertion of \<p 2 . The (3 functions are computed from 



of \ du 



f3 v (u,v) = pL 



dfj, 



i 

UO,VO 



(A2) 



u ,v 

while the RG functions r/^ and r\ t associated with the critical exponents are obtained from 

9 log Z^t 



(A3) 

UO,VQ 



d\og n 

The /5-functions have a simple dependence on d, indeed 

(3 U = (d-4)u + B u (u,v), (3 v = (d-4)v + B v (u,v), (A4) 

where the functions B u (u,v) and B v (u,v) are independent of d. Also the RG functions rj^t 
are independent of d. The standard critical exponents are related to rj^t by 

77 = t^(m*, v*), v=[2 + Vt (u*, v*) - V(/> (u*, v*)]" 1 ■ (A5) 

We report the five-loop series [17] for the cases N = 2 and 3. The series for N = 2 are: 

B u (u, v) = 2 u 2 - \uv + \v 2 - f m 3 + ^u 2 v - f A uv 2 + ^v 3 + 6.95758 m 4 - 2.95970 u 3 v 
+3.49036 u 2 v 2 - 1.45080 uv 3 + 0.0939364 v 4 - 33.3119 m 5 + 18.9022 u A v - 25.3312 u 3 v 2 
+14.5844 u 2 v 3 - 3.31596 uv 4 + 0.260717*; 5 + 197.427m 6 - 140.525 u 5 v + 211.453 u 4 v 2 
-152.787m 3 v 3 + 56.2903 uV - 10.8790mw 5 + 0.767937 v 6 , 

B v (u,v) = 2uv - \v 2 - fu 2 v + §uv 2 - £v 3 + 10.0721 u 3 v - 8.38186 u 2 v 2 + 2.73239 mv 3 
-0.272799 - 53.1466 u 4 v + 56.7468 m 3 ^ 2 - 29.2643 u 2 v 3 + 7.15431 m 4 - 0.598922 v 5 
+341.414 u 5 v - 444.234 mV + 311.112 u 3 v 3 - 120.254 mV + 22.4575mw 5 - 1.52571 v 6 , 

vA u , v ) = n« 2 - h uv + i8 y2 ~ + h u2y - m uv2 + ih v3 

+0.112847 m 4 - 0.0752315m 3 v + 0.0998264 mV - 0.041956m*; 3 + 0.00126591 v 4 - 0.41016m 5 
+0.3418 m\ - 0.466487 uV + 0.271787m V - 0.0630416 uv 4 + 0.00479535 v 5 , 

Vt (u, v) = -u + \v + \u 2 - \uv + \v 2 - |m 3 + §u 2 v - f^uv 2 + ^v 3 
+4.13719M 4 - 2.75813m 3 w + 2.95416mV - 1.18534mv 3 + 0.121022w 4 - 18.2814m 5 
+ 15.2345 u 4 v - 19.0871 u 3 v 2 + 10.7601 mV - 2.69595 uv 4 + 0.244579 v 5 . 

The series for N = 3 are: 

B u (u, v) = \u 2 - \uv + \v 2 - |m 3 + fu 2 v - §uv 2 + ^v 3 + 9.07446 m 4 - 6.28514 m 3 v 
+7.35688 m 2 v 2 - 2.98029 u v 3 + 0.149100 v 4 - 46.7683 m 5 + 42.8264 u 4 v - 56.8289 u 3 v 2 
+32.0271 mV - 6.92247 uv 4 + 0.503549 v 5 + 296.166 m 6 - 338.69 u 5 v + 506.708 mV 
-362.366 mV + 131.029 mV - 24.7085 uv 5 + 1.71145 v 6 , 
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B v (u,v) = 2uv - \v 2 - f u 2 v + f uv 2 - ±v 3 + 11. 1573m 3 w - 8.36797m 2 v 2 + 2.43934mv 3 
-0.264816 v 4 - 62.5357 m 4 v + 62.5357 uV - 30.8196 uV + 7.59503 mw 4 - 0.660444 
+422.235 m 5 v - 527.794 mV + 364.748 mV - 141.636 mV + 25.7555 mv 5 - 1.63935 v 6 , 

^(«, u) = \u 2 - ±uv + ±v 2 - j^U 3 + j^u 2 v - ^uv 2 + j^V 3 

+0.165895 m 4 - 0.165895 u 3 v + 0.228588 «V - 0.0935571 uv 3 - 0.00108507 v 4 - 0.709208m 5 
+0.88651 u 4 v - 1.19506 mV + 0.674671 mV - 0. 144646 uv 4 + 0.00961789 v 5 , 

r] t (u, v) = -§m + \v + §m 2 - \uv + - §m 3 + fu 2 v - luv 2 + |f^ 3 
+7.00217m 4 - 7.00217m 3 w + 7.32279 uV - 2.78612mw 3 + 0.227593w 4 - 33.6900m 5 
+42.1125m 4 w - 51.8428 uV + 28.3540 uV - 6.77118mv 4 + 0.599424w 5 . 

The critical exponents associated with the chiral degrees of freedom can be determined 
from the RG dimension of the chiral operator C c kdi{x) = 4>ck{x)(/)di{x) — </)ci{x)4>dk{x)- For 
this purpose, we computed the renormalization function Z c (u,v) by requiring Z C T^ C ' 2 ^ to be 
finite when expressed in terms of u and v. Here r^' 2 ) is the one-particle irreducible two-point 
function T^' 2 ) with an insertion of the operator C c kdi- Then, one defines the RG function 



d log Z c 

Vc{u,v) 



d log a , 

The resulting series for = 2 and N = 3 cases are respectively 

Vc (u, v) = -§« + \v + ^u 2 - ^uv + ±v 2 - |§m 3 + %u 2 v - %uv 2 + jgv 3 
+1.51005 u 4 - 3.99379 u 3 v + 3.36615 mV - 1.05138 uv 3 + 0.0947443 v 4 - 6.32935 m £ 
+ 19.2654 u 4 v - 21.2447 mV + 10.7984 mV - 2.47513mw 4 + 0.197296w 5 , 

and 



(A6) 



rj c (u, v) = -§m + I v + \u 2 - luv + > 2 - §;u 3 + §m 2 v - ffW + ^v 3 
+1.93436 m 4 - 5.12306 m 3 v + 4.24956 mV - 1.26059 uv 3 + 0.105233 w 4 - 8.73201m 5 
+26.6393 u 4 v - 29.198 mV + 14.5119 mV - 3.18347mv 4 + 0.246754 v 5 . 

The chiral crossover exponent C can be determined by using the RG scaling relation 

C = v [2 + 77 c (m*, u*) - ^(m*, v*)] . (A7) 

We have computed r/ c (u, v) for generic values of N. Thus, we have been able to compute the 
expansion of C in powers of e for any N. We have compared the result with the large-A?" 
expression of Ref. [82], finding full agreement. 

APPENDIX B: MEDIUM-RANGE MODELS 

In this Appendix we consider a rf-dimensional theory on a hypercubic lattice with general 
Hamiltonian 
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^ = -§E J p(^)E + E ^*)> ( B1 ) 

where <£>®| are M x N matrices, V(<^ x ) is an 0(M)(g)0(iV) invariant function, and the 
sums over x and y are extended over all lattice points. The coupling J p (x) depends on a 
parameter p. For instance, one may take the explicit form (3.9), but this is not necessary for 
the discussion that will be presented below. Indeed, one can consider more general families 
of couplings, as discussed in Sec. 3 of Ref. [63]. The relevant property is that J p (x) couples 
all lattice points for p — > oo, i.e., that for p — > oo one recovers a mean-field theory. 

Models like (Bl) are called medium-range models and admit an interesting scaling limit 
called critical crossover limit. If R parametrizes the range of the interactions, cf. Eq. (3.10), 
and f3 c {R) is the critical temperature as a function of R, then for R — > oo, t = ((3 C (R) — 
P)/Pc(R) — > at fixed t = R 2d ^ 4 ~ d H, critical quantities show a scaling behavior. For 
instance, the susceptibility x{fii R) an d the correlation length £((3,R) scale [62] according 
to Eq. (3.11). The functions f x (t) and fe(t) are directly related to the crossover functions 
F x (t,s) and F^(t,s) computed in field theory, cf. Eq. (3.12) [63,64]. The purpose of this 
Section is the computation of the nonuniversal constants p x , //£, a, and s. Interestingly 
enough, if the range R is defined according to Eq. (3.10), they do not depend on the explicit 
form of the coupling J(x), but only on the potential V(ip). The dependence on J(x) is 
effectively encoded in the variable R. 

The calculation can be done by a straightforward generalization of the results of Ref. [63]. 
Following Sec. 4.1 of Ref. [63] we first perform a transformation [83,84] — we use matrix 
notation and drop the subscript p from J p {x) to simplify the notation — rewriting [85] 



fa) = (detW)-™/* J exp (-fa* + , (B2) 



where <p is an iV x M matrix. Then, we define a function A(<f>) by requiring 



ze 



where z is a normalization factor ensuring A(0) = 0. We need to compute the small-0 
expansion of A(<f>). For this purpose we define the integrals 

W = / d NM <pe- y <»\<fT, (B4) 

» N M 

hn,2 = J d NM V e- V ^ £ £ W, (B5) 

ab cd 

and / n>m = I n ,m/Io,i- Then, a straightforward calculation gives 

N M 



A^ = a -f^+ a -f^r+ a -f 



^ Yl ( t )aC( t )ad( t )bC( t )bd ) - 



i2\2 



ab cd 



+ O(0 6 ), (B6) 



where 
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52,1 = ^Ai, (B7) 

- = HUa + 2/4,2) ^_f 2 ms^ 

4,1 MiV(M + 2)(iV + 2) N 2 M 2 ' 1 ' 1 ; 

- (MiV + 2)/ 4 , 2 -(M + iV + l)/ 4 , 1 

4 ' 2 MN(M -1)(N -1)(M + 2)(N + 2)' { } 

The expansion of A((f>) is all we need to compute the critical crossover limit. Indeed, it is 
possible to show that if correlations are directly related to correlations [63]. For instance, 

■ v»> = ~ (^) + (•?-')„ <*. ■ *•>■ (bio) 

WZ 

In the critical crossover limit, the first term in the right-hand side represents a subleading 
correction and can be neglected. This equation implies that x<p = ^2 \X<t>i where we have 
used the fact that P C (J2 X J ( x )) = V s 2,i + 0(R~ d ). 

Moreover, <f> correlations can be computed by using the following continuum theory 



N M 

<P ac <P ad <P bc <P bd ) - (</>■ 



< ^(a ) 2 \ 2 _ ^1 
4! 4! 



N M 

2\2 



ab cd 



(Bll) 



with a proper mass renormalization that is discussed in detail in Ref. [63]. This identification 
gives necessary conditions in order to observe a second-order 0(M)<g>0(iV) transition. The 
bare parameters should belong to the stability region, which implies a 4) i < and a 41 — 
\a^2 < 0, and should be in the attraction domain of the 0(M)®0(N) FP (assuming it 
exists), which implies (at least) a 4) 2 < 0. 

The above-reported results allow us to compute the nonuniversal constants appearing in 
Eq. (3.12) [86]. We obtain 

a 5 

= «2,i/4 = ( B12 ) 

a 4,l 

a = 5^, (B13) 
<i 

s = ^. (B14) 
a 4j i 

In Sec. IV we will be interested in the specific potential (3.7) for N = M = 2. In this case 



'2n,l 



vr 2 y^dy £ dx exp (-ry - §y 2 + ^Vy 2 ) , (B15) 



while I 2n ,2 can be obtained by taking derivatives of I 2n ,i with respect to V and U . Exact 
results can be obtained for V , U — > at fixed s L = V /?7o- It is easy to verify that 
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s m f(sL)/Uo — > if r < and s — > sl if r > 0, as expected in a weakly coupled system. 
We can also compute the nonuniversal constants and, in particular, s for the parameter 
values used in the MC simulation. For instance, we can consider two cases: (a) A 22 = 2 
that corresponds to r = —1, U = 36, and V = 24 (cf. Eq. 4.3); (b) A 22 = 7/5, i.e. 
r = -1, U Q = 144/5, and V = 48/5. For case (a) we obtain a 4 ,i = -0.0370438 and 
a«4 5 2 = —0.0336483, that satisfy the necessary conditions reported above. Correspondingly, 
s = 0.908337, a = 2.51497, fi x = 0.609569, and ^ = 3.22122. For case (b) we obtain 
a 4>1 = -0.0412887, a 4 , 2 = -0.0175521, s = 0.425106, a = 2.82353, fi x = 0.743715, and 
/j£ = 3.27408. It is interesting to note that for the family of Hamiltonians considered in 
Sec. IV — those with A 4 = 1 — we always have s > s L) the difference s — s L decreasing as 
sl — > 0. The exact mapping between the lattice model and the FT one also allows us to 
determine bounds on A* 22 for A4 — 1. Since s = 1 corresponds to A22 ~ 2.13 (sl ~ 0.72) 
and s = 2 to A 22 m 4.17 (s L w 1.23), the FT bound 1 < s* < 2 implies 2.13 < A* 2 < 4.17 
or, equivalently, 0.72 < s* L < 1.23. 
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